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Abstract 

We expand iterative numerically-exact influence functional path-integral tools and present a 
method capable of following the nonequilibrium time evolution of subsystems coupled to multi- 
ple bosonic and fermionic reservoirs simultaneously. Using this method, we study the real-time 
dynamics of charge transfer and vibrational mode excitation in an electron conducting molecular 
junction. We focus on nonequilibrium vibrational effects, particularly, the development of vibra- 
tional instability in a current-rectifying junction. Our simulations are performed by assuming large 
molecular vibrational anharmonicity (or low temperature). This allows us to truncate the molec- 
ular vibrational mode to include only a two-state system. Exact numerical results are compared 
to perturbative Master equation calculations demonstrating an excellent agreement in the weak 
electron-phonon coupling regime. Significant deviations take place only at strong coupling. Our 
simulations allow us to quantify the contribution of different transport mechanisms, coherent dy- 
namics and inelastic transport, in the overall charge current. This is done by studying two model 
variants: The first admits inelastic electron transmission only, while the second one allows for both 
coherent and incoherent pathways. 
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FIG. 1. Left Panel: Generic setup considered in this work, including a subsystem (S) coupled to 
multiple fermionic (F) and bosonic (B) reservoirs. Right panel: Molecular electronic realization 
with two metals, L and R, connected by two electronic levels, D and A. Electronic transitions in this 
junction are coupled to excitation/de-excitation processes of a particular, anharmonic, vibrational 
mode that plays the role of the "subsystem". This mode may dissipate its excess energy to a 
secondary phonon bath B. 



I. INTRODUCTION 



Following the quantum dynamics of an open-dissipative many-body system with multi- 
ple bosonic and fermionic reservoirs in a nonequilibrium state, beyond the linear response 
regime, is a significant theoretical and computational challenge. In the realm of molecular 
conducting junctions, we should describe the out-of-equilibrium dynamics of the molecu- 
lar unit while handling both electrons and molecular vibrations, accounting for many-body 
effects such as electron-electron, phonon-phonon and electron-phonon interactions. Given 
this complexity, studies in this field are mostly focused on steady-state properties, using 
e.g., scattering theory Q-Q, while ignoring vibrational nonequilibrium effects. Perturbative 
treatments (in either the molecule-leads coupling parameter or the electron-phonon interac- 
tion energy) are commonly used, includingthe nonequilibrium Green's function technique 
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4H8| and Master equation approaches [5|, |9 



13| . For following the real-time dynamics of 



such systems, involved methods have been recently developed, e.g., semiclassical approaches 

HQ. 

In this work, we extend numerically-exact path-integral methods, and follow the dynamics 
of a subsystem coupled to multiple out-of-equilibrium bosonic and fermionic reservoirs. The 
technique is then applied on a molecular junction realization, with the motivation to address 
basic problems in the field of molecular electronics. Particularly, in this work we consider the 



dynamics and steady-state properties of a conducting molecular junction acting as a charge 
rectifier. A scheme of the generic setup and a particular molecular junction realization are 
depicted in Fig. [TJ 

The time evolution scheme developed in this paper treats both bosonic and fermionic 
reservoirs. This is achieved by combining two related iterative path-integral methods: (i) 
The quasi-adiabatic path-integral approach (QUAPI) of Makri et al. [16], applicable for the 
study of subsystem-boson models, and (i) the recently developed influence-functional path- 



integral (INFPI) technique 17| . able to produce the dynamics of subsystems in contact 
with multiple fermi baths. The latter method (INFPI) essentially generalizes QUAPI. It 
relies on the observation that in out-of-equilibrium (and/or finite temperature) situations 
bath correlations have a finite range, allowing for their truncation beyond a memory time 
dictated by the voltage-bias and the temperature. Taking advantage of this fact, an iterative- 
deterministic time-evolution scheme can be developed, where convergence with respect to 
the memory length can in principle be reached. 

The principles of the INFPI approach have been detailed in [l^], where it has been 
adopted for investigating dissipation effects in the nonequilibrium spin-fermion model and 
charge occupation dynamics in correlated quantum dots. Recently, it was further utilized 
for examining the effect of a magnetic flux on the intrinsic coherence dynamics in a double 



quantum dot system 181 ] . and for studying relaxation and equilibration dynamics in finite 



metal grains 19 ]. 



Numerically-exact methodologies are typically limited to simple models; analytic results 
are further restricted to specific parameters. The Anderson-Holstein (AH) model has been 
studied extensively in this context. In this model the electronic structure of the molecule 
is represented by a single spinless electronic level, with electron occupation on the dot 
coupled to the displacement of a single oscillator mode, representing an internal vibration. 
This vibration may connect with a secondary phonon bath, representing a larger phononic 
environment (internal modes, solvent). The AH model has been simulated exactly with the 



secondary phonon bath, using a a real-time path-integral Monte Carlo approach 20j, and 



by extending the multilayer mu 



ticonfiguration time-dependent Hartree method to include 



fermionic degrees of freedom 21]. More recently, the model has been simulated by adopting 



the iterative-summation of path-integral approach 
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In this paper, we examine a variant of the AH model, the Donor (D)-Acceptor (A) elec- 



tronic rectifier model 



251 ] . This model incorporates nonlocal electron- vibration interactions: 



electronic transitions between the two molecular states, A and D, are coupled to a particular 
internal molecular vibrational mode. Within this simple system, we are concerned with the 
development of vibrational instability: Significant molecular heating can take place once 
the D level is lifted above the A level, as the excess electronic energy is used to excite the 
vibrational mode. This process may ultimately lead to junction instability and breakdown 



26l |. We have recently studied a variant of this model (excluding direct D-A tunneling ele- 



ment), using a Master equation method, by working in the weak electron-phonon coupling 
limit. {27]]. An important observation in that work has been that since the development of 
this type of instability is directly linked to the breakdown of the detailed balance relation 
above a certain bias (resulting in an enhanced vibrational excitation rate constant, over 
relaxation), it suffices to describe the vibrational mode as a truncated two- level system. In 
this picture, population inversion in the two-state system evinces on the development of 
vibrational instability. 



Our objectives here are threefold: (i) To present a numerically-exact iterative scheme 
for following the dynamics of a quantum system driven to a nonequilibrium steady-state 
due to its coupling to multiple bosonic and fermionic reservoirs, (ii) To demonstrate the 
applicability of the method in the field of molecular electronics. Particularly, to explore the 
development of vibrational instability in conducting molecules, (iii) To evaluate the per- 
formance and accuracy of standard-perturbative Master equation treatments, by comparing 
their predictions to exact results. Since Master equation techniques are extensively used 
for explaining charge transfer phenomenology, scrutinizing their validity and accuracy is an 
important task. 



The plan of the paper is as follows. In Sec. [IT] we introduce the path-integral formalism. 
We describe the iterative time evolution scheme in Sec. IHH by exemplifying it to the case of 
a spin subsystem. Sec. IIVI describes a molecular electronics application, and we follow both 
electrons and vibrational dynamics in a dissipative molecular rectifier. Sec. |V] concludes. 
For simplicity, we use the conventions H = 1, electron charge e = 1, and Boltzmann constant 



II. PATH-INTEGRAL FORMULATION 



We consider a multi- level subsystem, with the Hamiltonian Hs, coupled to multiple 
bosonic (B) and fermionic (F) reservoirs that are prepared in an out-of-equilibrium initial 
state. The total Hamiltonian H is written as 

H = H s + H B + H F + V SB + V SF . (1) 

In the energy representation of the isolated subsystem, its Hamiltonian can be written as 

#s = ^e s |s>(s| + ^^| S >( S / |. (2) 

S Sy^S' 

The Hamiltonian Hp may comprise of multiple fermionic baths, and similarly, Hb may 
contain more than a single bosonic reservoir. The terms Vsf and Vsb include the coupling 
of the subsystem to the fermionic and bosonic environments, respectively. Coupling terms 
which directly link the subsystem to both bosonic and fermionic degrees of freedom are not 
included. However, Vsb and Vsf may contain non- additive contributions with their own 
set of reservoirs. For example, Vsf may admit subsystem assisted tunneling terms, between 
separate fermionic baths (metals), see Fig. [TJ 

We are interested in the time evolution of the reduced density matrix ps(t). This quantity 
is obtained by tracing the total density matrix p over the bosonic and fermionic reservoirs' 
degrees of freedom 

ps(t) = Tr B Tr F [ e - im p(0)e iHt ] . (3) 

We also study the dynamics of certain expectation values, for example, charge current and 
energy current. The time evolution of an operator A can be calculated using the relation 

(A(t)} = Tt\p(p)A(t)} 

= lim ^-Tr \p(0)e lHt e XA e- tHt ] . (4) 

Here A is a real number, taken to vanish at the end of the calculation. When unspecified, the 
trace is performed over the subsystem states and all the environmental degrees of freedom. In 
what follows, we detail the path-integral approach for the calculation of the reduced density 
matrix. Section IIIIEI presents expressions useful for time-evolving expectation values of 
operators. 
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As in standard path-integral approaches, we decompose the time evolution operator into 
a product of N exponentials, e lHt = (e lH&t ^ N where t = N5t, and define the discrete time 
evolution operator Q = e lH5t . Using the Trotter decomposition, we approximate Q by 

G ~ QtQbQsQbQt-, (5) 

where we define 

Q T = e i( H F+V S F)5t/2 = e i(H B +V SB )St/2 

Qs = e iHs5t (6) 

Note that the breakup of the subsystem-bath term, e i(HF+v SF +H B +H SB )6t/2 ^ g B g T ^ j s exac t 
if the commutator \Vsb, Vsf] vanishes. This fact allows for an exact separation between the 
bosonic and fermionic influence functionals, as we explain below. This commutator nullifies 
if the fermionic and bosonic baths couple to commuting subsystem degrees of freedom, for 
example, Vsb oc \s)(s\ and Vsf oc \s')(s'\. 

As an initial condition, we assume that at time t = the subsystem and the baths are 
decoupled, p(0) = ps(0) ® ps <S> Pf, and the baths are prepared in a nonequilibrium (biased) 
state. For example, we may include in Hp two Fermi seas that are prepared each in a 
grand canonical state with different chemical potentials and temperatures. The overall time 
evolution can be represented by a path-integral over the subsystem states, 



(s%\Ps{t)\s N ) 

= - Yl Tl B^F (sN\G j \sN-l)( S N-l\G j \ S N-2)--(4 l S ) •••( S iV-2l^l S iv-l) ( S JV-1 l^l S Jv) 



°0 "1 "N-l 



(7) 

Here sf represents the discrete path on the forward (+) and backward (— ) contour. The 
calculation of each discrete term is done by introducing four additional summations, e.g., 

(8) 

We substitute Eq. ([8]) into Eq. ([7]), further utilizing the factorized subsystem-reservoirs 
initial condition as mentioned above, and find that the function under the sum can be 
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written product of separate terms, 



= E E E E E J *( m± > n± ' 4 M**, f ± , g ± )^(f ± , m ± , n^, g 

f=t 



(9) 



Here Is follows the subsystem (Hs) free evolution. The term Ip is referred to as a fermionic 
"influence functional" (IF), and it contains the effect of the fermionic degrees of freedom on 
the subsystem dynamics. Similarly, the bosonic IF, describes how the bosonic degrees 
of freedom affect the subsystem. Bold letters correspond to a path, for example, m 1 * 1 = 
{/71q , mf, mf ! _ 1 }. We also define the path = {sq , sf, and the associate 
path which covers N + 1 points, s /=t = {sf , sf , s%_ 1: sf}. Given the product structure of 
Eq. (Q, the subsystem, bosonic and the fermionic terms can be independently evaluated, 
while coordinating their path. Explicitly, the elements in Eq. (Q are given by 



Is 
If 



\4\Ps{fy\s )nk=o,...,N-i(m k \Gs\n k )(n+\gl\m+) 



s n I @f I 9n- i ) (/iv-i I I S N- 1 ) • • • 



1 1„+ 



x (4\GU9o){fo\GF\4 )Pf(s \g F \f )(9 \Qf\s 1 }... 



X (s 



N- 



-l\GF\fN-x)(9N-l\0F\ 8 N) 
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Tr 



B 



(^_il^sK_ 1 )(m^_ 1 |^|/^_ 1 )... 
x (^|£i?K>( m o \G B \fo)PB(f \GB\m )(n Q \Gb\9o)- 



(10) 



The dynamics in Eq. (Q can be retrieved by following an iterative scheme, by using the 
principles of the INFPI approach [l7i]. In the next section we illustrate this evolution with 
a spin subsystem. 



III. ITERATIVE TIME EVOLUTION SCHEME 



We consider here the spin-boson-fermion model. It includes a two-state subsystem that 
is coupled through its polarization to bosonic and fermionic reservoirs. With this relatively 
simple model, we exemplify the iterative propagation technique, see Sees. IIII AlllTTEl Rel- 
evant expressions for a multi-level subsystem and general interaction form are included in 
Sec. IfflFl 
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A. spin-boson-fermion model 



The spin-fermion model, with a qubit, spin, coupled to a fermionic bath is kindred to 
the eminent spin-boson model, describing a qubit interacting with bosonic environment. It 



is also related to the Kondo model 



only lacking direct coupling of the reservoir degrees 



of freedom to spin-flip processes. It provides a minimal setting for the study of dissipation 



and decoherence effects in the presence of nonequilibrium reservoirs 



29|-|32|. Here we put 



together the spin-boson and the spin-fermion models, and present it in the general form, 

H s = Aa x + Ba z , 

H f = Y1 e A c i + £ v lA c i' 

hi' 

p Pip 1 

v SB = £ ( b l + b p) + a * E OJv- (ii) 

p p,p' 

The subsystem includes only two states, with an energy gap 2B and a tunneling splitting 
2A. This minimal subsystem is coupled here through its polarization to a set of boson 
and fermion degrees of freedom, where a z and a x denote the z and x Pauli matrices for 
a two-state subsystem, respectively. b p stands for a bosonic operator, to destroy a mode 
of frequency u p , similarly, Cj is a fermionic operator, to annihilate an electron of energy tj 
(we assume later a linear dispersion relation). In this model, spin polarization couples to 
harmonic displacements, to scattering events between electronic states in the metals (fermi 
reservoirs), and to scattering evens between different modes in the harmonic bath. Since the 
commutator between the interaction terms vanish, [Vsf, Vsb] = 0, the separation between 
the bosonic and fermionic IFs is exact. Moreover, since the fermionic and bosonic operators 
couple both to a z , we immediately note that f k — s k , nrj: = f k , = g k and g k = s k+1 . 
Eq. OU) then simplifies to 

(s + N \ps(t)\s- N ) = ^/ s (s' ± )/ F (s /± )/ B (s' ± ), (12) 

s± 

where we recall the definitions of the paths = {sq , sf, s^-_ 1 } and s /=t = {sq , sf, s^_ 1; 
The subsystem evolution and the IFs are now given by 



Is(s f± ) 
I F (s'±) 
where 



(s+\p s (0)\s )K(s%, st)K(sf, 4) 

iW B (s+)8t/2 -iWsis+^St -iW B {s+)St/2 iW B (So)5t/2 (n iW B (s^_ 1 )5t JW B (s^)6t/2 



p -iW F (s+)St/2 p-iWpis+^St p -iW F {s+)St/2 pF ^W F {s^)8t/2 p iW F {s-_{)St p m F {s-)St/2 



(13) 



\ s k+i\ e 



-iH s St\ c + 



S k)( S k\ e 



iH s 6t\ 



'fe+l 



) 



(14) 



is the propagator matrix for the subsystem. We have also used the short notation W for 
bath operators that are evaluated along the path, 



W F (s) = H F +(s\V SF \s), 
W B (s) = H B + (s\V SB \s). 



(15) 



In the next sections we explain how we compute the bosonic and fermionic IFs. The former 
has a closed analytic form in certain situations. The latter is computed only numerically. 

B. Bosonic IF 



We present the structure of the bosonic IF in two separate models, corresponding to 
different types of subsystem-boson bath interactions. In both cases the bosonic bath is 
prepared in a canonical state of inverse temperature (3 ph = l/T ph , 

p B = e ~PvhH B /Tr B [e-^ hHB }. 



(16) 



Displacement interaction model, vL 3 , = and Cf , = 0. Given the remaining linear 
displacement-polarization interaction, an analytic form for the bosonic IF can be written, 
the so-called "Feynman- Vernon" influence functional (FV IF) [33] . In its time-discrete form, 
the bosonic IF is given by an exponent with pairwise interactions along the path 

N k 

Ib{sq,...,s^) = exp ~^2^2(st ~ s^)(Vk,k'S^ ~ V*k,k' s k') ■ ( 17 ) 

. k=0 k'=0 

The coefficients rj^k' are additive in the number of thermal baths, and they depend on these 
baths' spectral functions and initial temperatures 
are included in Appendix A. 



161 ] . For completeness, these coefficients 



Boson scattering model, £f = 0. The bosonic IF can now be computed numerically, by 



using the trace formula for bosons 



34] 



Tr B [e Ml e M2 ...e Mk ] = det[l - e mi e m2 ...e mh )- 1 . (18) 

Here is a single particle operator corresponding to a quadratic bosonic operator Mf. = 
J2p P >( m k)p,p'blb p >. Application of the trace formula to the bosonic IF (TTBT) leads to 

I B = Tr B [e^e M \..e Mk p B } 

= det{[/B + f B ] - e m ie m *...e m *f B y\ (19) 

The matrix I B is an identity matrix, and the function f B stands for the Bose-Einstein 
distribution, f B = \eP vhU — l] -1 . The determinant in Eq. ( |T9l can be evaluated numerically 
by taking into account L B modes for the boson bath. This discretization implies a numerical 
error. Generalizations, to include more that one bosonic baths, are immediate. 

C. Fermionic IF 

The fermionic IF is computed numerically since an exact analytic form is not known in 



the general strong coupling limit [29|-|31|. It is calculated by using the trace formula for 



fermions 



Tr F [e Ml e M2 ...e Mk ] = det[l + e mi e m2 ...e mfc ]. (20) 



Here is a single particle operator corresponding to a quadratic operator Mt = J2i j( m k)i,j c l c j- 
In the next section we consider a model with two Fermi seas, Hp = Hp + Hr, prepared in 
a factorized state of distinct grand canonical states, p F = Pl® Pr, with 

p v = e-M^-^A^/Tr^e-M",-/^,)^ U = L,R (21) 

Here @ v = 1/T U stands for an inverse temperature, and p v denotes the chemical potential of 
the v bath. Application of the trace formula to the fermionic IF in Eq. f lT5|) leads to 

I F = Tr F [e Ml e M2 ...e Mfc P F] 

= det{[/ L - f L \ ® [Ir - f R ] + e mi e m2 ...e m *[/L ® f R ]j. (22) 

The matrices I v are the identity matrices for the v = L,R space. The functions and Jr 
are the bands electrons' energy distribution, f v = [e^"( e_ ^) + l] -1 . The determinant in Eq. 
( |22|) can be evaluated numerically by taking into account L s electronic states for each metal. 
This discretization implies a numerical error. 
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D. The iterative scheme 



The dynamics described by Equation ( 1121) includes long-range interactions along the 
path, limiting brute force direct numerical simulations to very short times. The iterative 
scheme, developed in Ref. is based on the observation that in standard jionequilibrium 



situations and at finite temperatures bath correlations exponentially die [24J, |3l| , thus the IF 
can be truncated beyond a memory time r c = N s 5t, corresponding to the time where bath 
correlations sustain. Here N s is an integer, St is the discretized time step, and the correlation 
time r c is dictated by the bias and temperature. Roughly, for a system under a potential 
bias A/i and a temperature T, r c ~ max{l/T, 1/A/j,} [171 ] . By recursively breaking the 
IF to include terms only within r c , we reach the following (non-unique) structure for the 
a = B,F influence functional, 

J-aV^O J *1 J *2 J ••■> °N) 1 a\ b Q j 6 1 5 •••) b N s ) 1 a \ b l ) *2 ) "•> i, V s +l/-'o 1*2 > 6 3 ) •••) *iV s +2/-" 
X -^Q! S \ S N-N s i S N-N 3 +n "•' S Ar)> (23) 

where we identify the "truncated IF" , 14 , as the ratio between two IFs, with the numerator 
calculated with an additional time step, 

T (Ns)( Q b \ _ ^( s fc; s fc+n •••; s fc+JV a ) 

i «l s fc ' s fe+l' "■' s fc+iV i ,-l/ 

(24) 

The truncated IF is the central object in our calculations. For fermions, its numerator and 
denominator are separately computed using Eq. (l2~2~j) . The bosonic IF is similarly computed 
with the help of Eq. (119p when ^ = 0. In the complementary case, (p p , = and , = 0, 
the truncated-bosonic IF has a closed analytic form: Using Eq. ( ITTj) we find that it comprises 
only two-body interactions, of s^+Ns with the preceding spins, down to s&, 

k+N s 

( S fc+V a — S k+N s )( r lk+N s ,k'S^ / — V*k+N s ,k' s k' 



(Sfc, Sfc+i, s fc+J v s ) = exp 



(25) 



Based on the decompositions ( 124)) and ( 1251) . we time-evolve Eq. (fT2|) iteratively, by defining 
a multi-time reduced density matrix ps{sk, Sk+h •-■> s k+N s -i)- Its initial value is given by 

Ps(4, •••> 4j = J 5(4 , •-, s^ s )/ B (4, •-, sJ,)M s o> s%J. (26) 
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Its evolution is dictated by 



X ^F \ S k ! ) S k+N S )^B \ S k> S k+N a )- (2?) 

The time-local (t^ = k8t) reduced density matrix, describing the state of the subsystem at 
a certain time, is reached by summing over all intermediate states, 

Ps(t k )= Ps(4-n.+v-,4)- (28) 

s ± s ± 

The bosonic and fermionic IFs may be (and often this is the case) characterized by different 
memory time. Thus, in principle we could truncate the fermionic IF to include iVf terms, 
and the bosonic IF to include elements. However, the efficiency of the computation is 
dictated by the longest memory time, thus, for convenience, we truncate both IFs using the 
largest value, identified by N s . 

By construction, this iterative approach conserves the trace of the reduced density matrix, 



ensuring the stability of the iterative algorithm to long times [16(|. This property can be 
inferred from Eqs. ffl2|) and f fl3|) . by using the formal expressions for the truncated IFs, Eq. 
( I23|) and (f24|) . To prover this property, we trace over the reduced density matrix at time t, 
identifying s N = s% = s^, 

Trs\ps(t)] = Y,MPs(t)\8 N ) 

SN 
s'± 

Using the cyclic property of the trace, we note that both the fermionic and bosonic IFs 
are independent of sn, when = sjf. Therefore, the summation over the sn coordinate 
reduces to a simple sum which can be performed using the completeness relation for the 
subsystem states, resulting in 

^( S Ar|e-^ 5t |4-i)(^-il^|^) = 5(4-i - ^-i)- (29) 
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Iterating in this manner we conclude that 



T±s\ps(t)] = ^2(s N \ps(t)\s N ) 



= ^/ s (s /± )J J ,(s /± )/ fl (s /± )5(4 - s N )5(s 



.+ 

N-l 



s 1 )S(s, 



.+ 





= ^o|ps(0)|s > = Tr 5 [p 5 (0)] 



(30) 



We emphasize that the trace conservation is maintained even with the use of the truncated 
form for the IFs. Moreover, it holds irrespective of the details of the bath and the system- 
bath interaction form. It is also obeyed in the more general case, Eq. Q. Equation ( 127)) 
[and its generalized form, Eq. (I34p below], describe a linear map. Its fixed points are stable 
if the eigenvalues of the map have modulus less than one, which is the case here. Thus, our 
scheme is expected to approach a stationary-state in the long time limit. 

E. Expectation values for operators 

Besides the reduced density matrix, we can also acquire the time evolution of several 
expectation values. Adopting the Hamiltonian pip , we illustrate next how we achieve 
the charge current behavior. For simplicity, we consider the case with only two fermionic 
reservoirs, v = L,R. The current operator, e.g., at the L bath is defined as the time 
derivative of the number operator. The expectation value of this current is given by 



We consider the time evolution of the related exponential operator e L , with A a real 
number that is taken to vanish at the end of the calculation, 



As before, the initial condition is factorized at t — 0, p(0) = ps(0) <8> Pb <S> Pf- The trace is 
performed over subsystem and reservoirs degrees of freedom. By following the same steps 
as in Eqs. (J3])-(J7j), we reach the path-integral expression 




(31) 



(N L (t)) = TtipN L (t)} 




lim — Tr \p(0)e iHt e XNL e- lHt ] . 



(32) 
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5 A l b N~l 



x (4Ip(0)I s o)---( s jv-2I^I%-i)( s 7v-iI^I s ^)J- (33) 
Factorizing the time evolution operators using Eq. (jSJ), we accomplish the compact form 

{e XN L ( t)) = J2ls(s' ± )lB(s' ± )j F (s' ± )S(s + N - S N )- 
s'± 

The terms is and Is are given in Eq. ffT5|) . The fermionic IF accommodates an additional 
exponent, 

J F (s /=t ) = Tr F \ e ^L e -iW F {8+)st/2 e -iw F (8+_ 1 )6t^^ e -iw F {a+)st/2^ 

We can time evolve the operator (e XNL ) by using the iterative scheme of Sec. III.D, by 
truncating the bosonic and fermionic IFs up to the memory time r c = N s 5t, for several 
values of A. We then take the numerical derivative with respect to A and t, to attain the 
charge current itself. 

The approach explained here could be used to explore several fermionic operators, for 
example, the averaged current j av = (Jz — jo)/ 2. The minus sign in front of ]r originates 
from the sign notation, with the current defined positive when flowing L to R. The im- 
plementation of a heat current operator, describing the heat current flowing between two 
bosonic reservoirs, requires first the derivation of an analytic form for the bosonic IF, an 
expression analogous to the FV IF, and the subsequent time discretization of this IF, to 
reach an expression analogous to (|17[) . 



F. Expression for multilevel subsystems and general interactions 

So far we have detailed the iterative time evolution scheme for the spin-boson-fermion 
model (TTTT) . The procedure can be extended, to treat more complex cases. Based on the 
general principles outlined in Sec. IHI D\ one notes that the path-integral expression (Q can 
be evaluated iteratively by generalizing Eq. (127)) to the form 

^^Ps( v k 5 • ) v k+N s -l)-^( m k+N s ^ U k+N S )^F S \ S k J fk ' 9k ' S k+N a i fk+N s i 9k+N a ) 

X ^B \fki9k> m ki n k>--->fk%N a >9k+N 3 i m k+N s i n k+N a ) (34) 
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where we compact several variables, v k = {s k , f k ,g k , m k ,n k }. It should be noted that in 
cases when the IF is time invariant, as in the molecular electronics case discussed below, 
one needs to evaluate Ig and only once, then use the saved array to time-evolve the 
auxiliary density matrix. 
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FIG. 2. Molecular electronic rectifier setup. A biased donor-acceptor electronic junction is coupled 
to an anharmonic mode, represented by the two-state system with vibrational levels |0) and |1). 
This molecular vibrational mode may further relax its energy to a phononic thermal reservoir. This 
process is represented by a dashed arrow. Direct electron tunneling element between D and A is 
depicted by a dotted double arrow. Top: A/x > 0. In our construction both molecular electronic 
levels are placed within the bias window at large positive bias, resulting in a large (resonant) 
current. Bottom: At negative bias the energy of A is placed outside the bias window, thus the 
total charge current is small. 



IV. APPLICATION: MOLECULAR RECTIFIER 

The functionality and stability of electron-conducting molecular junctions are directly 
linked to hea ting and cooling effects experienced by molecular vibrational modes in biased 



situations 



35 



-411]. In particular, junction heating and breakdown may occur once the 
bias voltage exceeds typical molecular vibrational frequencies, when the electronic levels are 
situated within the bias window, if energy dissipation from the molecule to its environment 
is not efficient. 
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In this section we study the dynamics and steady-state behavior of electrons and a spe- 
cific vibrational mode in a molecular conducting junction serving as an electrical rectifier. 
The rectifier Hamiltonian is detailed in Sec. IIV Al In Sec. IIVBI we show that this model 
can be mapped onto the spin-boson- fermion Hamiltonian This allows us to employ the 
path-integral technique of Sec. II III for simulating the rectifier dynamics. The rectification 
mechanism is explained in Sec. IIV CI Relevant expressions of a (perturbative) Master equa- 
tion method are described in Sec. IIV to be compared to our path-integral based results 
in Sec. IIV El Convergence issues and computational aspects are discussed in Sec. IIVFI 



A. Rectifier Hamiltonian 

The D-A rectifier model includes a biased molecular electronic junction and a selected 
(generally anharmonic) internal vibrational mode which is coupled to an electronic tran- 
sition in the junction and to a secondary phonon bath, representing other molecular and 
environmental degrees of freedom. In the present study we model the anharmonic mode 
by a two-state system, and this model can already capture the essence of the vibrational 



instability effect [27|]. For a schematic representation, see Fig. [2j This model allows us to 
investigate the exchange of electronic energy with molecular vibrational heating, and the 
competition between elastic and inelastic transport mechanisms. Its close variant has been 



adopted in Refs. 



42 



2|-|44j for studying the thermopower and thermal transport of electrons 
in molecular junctions with electron-phonon interactions, within the linear response regime. 

We assume that the D molecular group is strongly attached to the neighboring L metal 
surface, and that this unit is overall characterized by the chemical potential //£,. Similarly, 
the A group is connected to the metal R, characterized by fiR. At time t = the D and 
A states are put into contact. Experimentally, the R metal may stand for an STM tip 
decorated by a molecular group. This tip is approaching the D site which is attached to the 
metal surface L. Once the D and A molecular groups are put into contact, electrons can flow 
across the junction in two parallel pathways: (i) through a direct D-A tunneling mechanism, 
and (ii) inelastically, assisted by a vibration: excess electron energy goes to excite the D-A 
vibrational motion, and vice versa. 

The rectifier (rec) Hamiltonian includes the electronic Hamiltonian H e i with decoupled 
D and A states, the vibrational, two-state subsystem H vib , electronic-vibrational coupling 
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Hi, a free phonon Hamiltonian E p h, and the coupling of this secondary phonon bath to the 
selected vibration, 

H r ec — H e i + H v ii, + Hj + Hph + Hyib-ph. (35) 

The electronic (fermionic) contribution E d attends for all fermionic terms besides the direct 
D and A tunneling term, which for convenience is included in Hj, 

H e i = H M + + Hx + E c 

H M = e d c\c d + e a c\c a 

H l = Y1 e ' c i c '! H °R = Y1 €rC r Cr - 

leL reR 

H C = ^Vi (c\c d + C^C^j + ^ V r (4 c a + 4 c r) • (36) 
I r 

Hm stands for the molecular electronic part including two electronic states, a donor D and 
an acceptor A. cj^ (cd/ a ) is a fermionic creation (annihilation) operator of an electron on 
the D or A sites, of energies e d>a . The two metals, v = L, R, are each composed of 
a collection of noninteracting electrons. The hybridization of the D state to the left (L) 
bath, and similarly, the coupling of the A site to the right (R) metal, are described by 
Hq- The vibrational Hamiltonian includes a special nuclear anharmonic vibrational mode 
of frequency ojq, 

H vi b = (37) 

The displacement of this mode from equilibrium is coupled to an electron transition in the 
system, with an energy cost k, resulting in heating and/or cooling effects, 

Ej = (na x + v da ) (c d c a + c\c d ^j . (38) 

Besides the electron-vibration coupling term, Ei further includes a direct electron tunneling 
element between the D and the A states, of strength v da . Electron transfer between the two 
metals can therefore proceed through two mechanisms: coherent tunneling and vibrational- 
assisted inelastic transport. 

The selected vibrational mode may couple to many other phonons, either internal to the 
molecules or external, grouped into a harmonic reservoir, 

p 

H vib - ph = cr x J2 ^ K + b p ) (39) 
p 
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The Hamiltonian H V Q,- P h corresponds to a displacement-displacement interaction type. 
The motivation behind the choice of the two-level system (TLS) mode is twofold. First, 



as we showed in Ref. 



271 ]. the development of vibrational instability in the D-A rectifier 



does not depend on the mode harmonicity, at least in the weak electron-phonon coupling 
limit. Since it is easier to simulate a truncated mode with our approach, rather than a 
harmonic mode, we settle on the TLS model. Second, while there are many studies where 



a perfectly harmonic mode is assumed, for example, see Refs. [20|, [21|, |24j , to the best of 
our knowledge our work is the first to explore electron conduction in the limit of strong 
vibrational anharmonicity. 



B. Mapping to the spin-boson- fermion model 

We diagonalize the electronic part of the Hamiltonian H e i to acquire, separately, the exact 
eigenstates for the L-half and R-half ends of H e i, 

H e i = H L + H R 

H L = ^eiajai, H R = y]e T ala r . (40) 

l r 

Assuming that the reservoirs are dense, their new operators are assigned energies that are 
the same as those before diagonalization. The D and A (new) energies are assumed to be 
placed within a band of continuous states,excluding the existence of bound states. The old 



operators are related to the new ones by 



45] 



Cd = ^ Xm, ci = r]i tV a v 

I V 



^■(2 — ^^^^^ hyO/ip j Cy — ^^^^^ ^~\t 7*' ^T -1 ' J ^4 1 ^ 



where the coefficients, e.g., for the L set, are given by 

Vl 



A, 



ei-e d - J2i> 



Vi,l> = $i,l> ~ ~ -—-7- ( 42 ) 

e« — €i> + id 

Similar expressions hold for the R set. It is easy to derive the following relation, 

E VJ -^ = PP J2— (43) 
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with the hybridization strength (vj is assumed real), 

T L (e) = 2nJ2^-ei). (44) 
l 

With the new operators, the Hamiltonian (135]) can be rewritten as 

H rec = ^2 e l a \ a l + 6 r a i a r + 
I r 

+ {Ka x + v da ) ^2 ^A r a}a r + \*\[alai 

l,r 

+ J2 u} p b l b p +a *Y,^( b l +b p)- ( 45 ) 

V V 

This Hamiltonian can be transformed into the spin-boson-fermion model of zero energy 
spacing, using the unitary transformation 

U jf a z U = a x , U^a x U = a Z} (46) 

with U = -j^ipx + &z)- The transformed Hamiltonian H rec = WH rec U includes a cr^-type 
electron-vibration coupling, 

Hrec = ^2 € l a \ a l + e r a l a r + ~^ a x 

I r 

+ {ko z + v da ) ^2 [A*A r a}a r + \* r \ia\ai 

l,r 

+ " P b ] p b p + a z J2 iv ( b l + h v) ■ ( 47 ) 

v v 

It describes a spin (TLS) coupled diagonally to two fermionic environments and to a single 
boson bath. One can immediately confirm that this Hamiltonian is accounted for by Eq. 
( ITTj) . To simplify our notation, we further identify the electronic- vibration effective coupling 
parameter 

e£ = k$k. (48) 

For later use we also define the spectral function of the secondary phonon bath as 

J ph (u) = nJ2& 2 ^-" P )- (49) 



v 



In our simulations below we adopt an ohmic function, 



J ph (u) = (50) 
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with the dimensionless Kondo parameter Kd, characterizing subsystem-bath coupling, and 
the cutoff frequency u c . 

As an initial condition for the reservoirs, we assume canonical distributions with the 
boson-phonon bath distribution following pg = e~P phHph /TiB[e~ l5phHph ] and the electronic- 
fermionic initial density matrix obeying pp = Pl®Pr-, with p u = e -Pv( H v-VvN u ) j r Yxp\e~^ v<yUl '^ iL ' 
v = L,R. This results in the expectation values of the exact eigenstates, 

(a\a v ) = 5i iV f L (ei), (a\a r ,) = 5 r>r ,f R (e r ), (51) 

where /z,(e) = [exp(/?i(e — pi)) + l] -1 denotes the Fermi distribution function. An analogous 
expression holds for The reservoirs temperatures are denoted by 1/(3 U ; the chemical 

potentials are p v . 
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FIG. 3. Left panel: Energy of the donor (full line) and acceptor states (dashed line). The dotted 
lines correspond to the chemical potentials at the left and right sides. Right panel: Damping rate 
K v tf,. The junction's parameters are T u = 1, (3 U = 200, k = 0.1, ujq = 0.2, and e^Ap, = 0) = —0.2, 
e a (Ap = 0) = 0.4. We used fermionic metals with a linear dispersion relations for the original H® 
baths and sharp cutoffs at ±1. All energy parameters are given in units of eV. 



C. Rectifying mechanism 

We now explain the operation principles of the molecular rectifier. In our construction the 
application of a bias voltage linearly shifts the energies of the molecular electronic levels, D 
and A. In equilibrium, we set e a < and > 0. Under positive bias, defined as pi — pr > 0, 
the energy of the acceptor level increases, and the donor level drops down, see Fig. |5J When 
both levels are buried within the bias window, the junction can support large currents. At 
negative bias the electronic level A is positioned above the bias window, resulting in small 
currents. For a scheme of the energy organization of the system, see Fig. [3] panel, left panel. 
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FIG. 4. Scheme of the vibrational mode excitation and relaxation processes. A full circle represents 
an electron transferred; a hollow circle depicts the hole that has been left behind. 

A generic mechanism leading to vibrational instabilities (and eventually junction rupture) 



in D-A molecular rectifiers has been discussed in Ref. [26[: At large positive bias, when the 
D state is positioned above the acceptor level, electron-hole pair excitations by the molecular 
vibration (TLS) dominate the mode dynamics. This can be schematically seen in Fig. HI 
The second-order perturbation theory rate constant, to excite the vibrational mode, while 
transferring an electron from L to R, k^Zti ■, overcomes other rates once the density of states 
at the left end is positioned above the density of states at the right side. This is the case at 
large positive bias, given our construction. The rate k^-^ is defined next, in Sec. IIVDI 



D. Master equation (vd a = 0) 

In the limit of weak electron-vibration coupling, once the direct tunneling term is ne- 
glected, t>da=0, it can be shown that the population of the truncated vibrational mode 
satisfies a kinetic equation (27]. 

Pi = - (kUo + k Uo) Pi + (fcJUi + k o^i) Po, 

Po+Pi = l- (52) 

The excitation (fco->i) an d relaxation (ki-+o) rate constants are given by a Fourier transform 
of bath correlation functions of the operators F e and F&, defined as 

F e = ^2(tf r a\a r + £r,l a l a l)> 

l,r 

F h = T,£( b l +b p)> ^ 
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to yield 



/oo 
e^-^ T Ti F [p F F e (r)F e (0)]dr 
-oo 

S € '- e ^ T Tr B \pBF h {T)F h {ti)\dT. 



A e — % a' 



OO 

oo 



(54) 



Here s — 0, 1 and e\ — e = Uq. The operators are given in the interaction representation, 



e.g., a\(t) = e iHLt a\e' iHLt . 



Phonon-bath induced rates. Expression (154"]) can be simplified, and the contribution of 
the phonon bath to the vibrational rates reduces to 



k 1^0 = r ph(^o)[/B(Wo) + 1], 
h,b — l.b p -^oP v h 



b 0^1 



(55) 



where = [e^ phUJ — 1] 1 denotes the Bose-Einstein distribution function. The damping 

rate is defined as T p h(uj) = 2J p h(uj), 



(56) 



For brevity, we ignore below the direct reference to frequency. 

Electronic-baths induced rates. The electronic rate constants (154)) include the following 



contributions 



L^R | uR^L 



satisfying 



k^ R = 27TK 2 \M 2 \K\ 2 h(ei)(l " f R (e r ))S(u; + e l - e r ) 

l,r 

k^ R = 27TK 2 M 2 |A r | 2 / L (Q)(l - f R (e r ))S(-uj + e t - e r ). 

hr 



(57) 



(58) 



Similar relations hold for the right-to-left going excitations. The energy in the Fermi function 
f u (e) is measured with respect to the (equilibrium) Fermi energy, placed at + [ir), and 
we assume that the bias is applied symmetrically, /i^ = — /ir. The rates can be expressed 
in terms of the fermionic v = L,R spectral density functions 



J u (e) = 2nKj2\\ 3 \ 2 5(e 3 -e) 



(59) 
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Using Eq. (142]) we resolve this as a Lorentzian function, centered around either the D or the 
A level, 



JAe) = k- 



e-e d ) 2 + r L (e) 2 /4 

r R (e) 



The electronic hybridization T u (e) is given in Eq. (jUJ). Using these definitions, we express 
the electronic rates [Eq. flo^]) ] by integrals (s, s'=0,l) 

1 r°° 

KZ = T~ / t 1 " ^'( e + ( s - s >°)] Ut)JAt + (« - O^o)de. (61) 

Observables. Within this simple kinetic approach, junction stability can be recognized 
by watching the TLS population in the steady-state limit: population inversion reflects on 



vibrational instability 



27|. Solving Eq. ( 152]) in the long time limit we find that 

h e 4- b b _|_ U e 4- h b '' 



px= ^ ^riir^a . Pd=i-Pi. (62) 



A related measure is the damping rate K vi b [26j, depicted in Fig. [3] panel (b). It is defined 



as the difference between relaxation and excitation rates, 

K V ib = fci-^o + ^i^o ~~ (^o-+i + ^o-+i) • (63) 

Positive K V ib indicates on a "normal" thermal-like behavior, when relaxation processes over- 
come excitations. In this case, the junction remains stable in the sense that the population 
of the ground state is larger than the population of the excited level. A negative value for 
K V ib evinces on the process of an uncontrolled heating of the molecular mode, eventually 
leading to vibrational instability and junction breakdown. 



In the steady-state limit, the charge current j, flowing from L to R, is given by 27] 



3 = Pi ~ k^ L ) + Po (k^ R - k^ L ) . (64) 

This relation holds even when the TLS is coupled to an additional boson bath. Note that in 
the long time limit the current that is evaluated at the left end ji is equal to j R . Therefore, 
we simple denote the current by j in that limit. 

Master equation calculations proceed as follows. We set the hybridization energy T v as 
an energy independent parameter, and evaluate the fermionic spectral functions J v {e) of Eq. 
(160]) . With this at hand, we integrate (numerically) Eq. fl6Tl) . and gain the fermionic-bath 
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induced rates. The phonon bath-induced rates (1551) are reached by setting the parameters of 
the spectral function J p h, to directly obtain T p h, see Eq. fl56|) . Using this set of parameters, 
we evaluate the levels occupation and the charge current directly in the steady-state limit. 
We can also time evolve the set of differential equations ( l5"2"j) . to obtain the trajectory pi t o(t). 




FIG. 5. Absolute value of the quantity Trp£f r . The figure was generated by discretizing the 
reservoirs, using bands extending from —D to D, D = 1, with Nl = 200 states per each band a 
linear dispersion relation and a constant density of states for the R reservoirs, with a constant 
density of states p = Nl/2D. Electron- vibration coupling is given by k = 0.1. 



E. Results 

We simulate the dynamics of the subsystem in the spin-boson-fermion Hamiltonian ( 14"T|) 
using the path-integral approach of Sec. Ill II In order to retrieve the vibrational mode 
occupation in the original basis in which Eq. (145]) is written, we rotate the reduced density 
matrix ps{t) back to the original basis by applying the transformation U = -j^ipx + °z)> 

Ps(t) = U Ps {t)U. (65) 

The diagonal elements of ps(t), correspond to the vibrational mode occupation, the ground 
state |0) and the excited state |1), 

Po(*) = (0|p s (t)|0) Pl (t) = (l\p s (t)\l). (66) 

As an initial condition we usually take ps(0) = |(— o x + / s ), J s is a 2 x 2 unit matrix. Under 
this choice, ps(0) has only its ground state populated. 
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Our simulations are performed with the following setup, displayed in the left panel of 
Fig. [3j In the absence of a bias voltage we assign the donor the energy = —0.2 and the 
acceptor the value e a = 0.4. These molecular electronic states are assumed to linearly follow 
the bias voltage. The right panel in Fig. [3] depicts the damping rate K v ^ in the absence of 
coupling to the phonon bath, as evaluated using the Master equation method. This measure 
becomes negative beyond A/x ~ 0.85, which corresponds to the situation where the (bias 
shifted) donor energy exceeds the acceptor by wo, Q — e a > wo; ojq = 0.2. This results in a 
significant exchange of electronic energy to heat, affecting junction's instability. 



1. Isolated mode 

We study the time evolution of the vibrational mode occupation using vja = (unless 
otherwise stated), further decoupling it from a secondary phonon bath, Kd=0. 

Electron-vibration interaction energy. The interaction energy of the subsystem (TLS) 
to the electronic degrees of freedom is encapsulated in the matrix elements £ ; F r = n\^\ r , 
see Eq. (j4"8j) . The strength of this interaction is measured by the dimensionless parameter 
Trp(ep)£f r , which connects to the phase shift experienced by Fermi sea electrons due to a 



scattering potential, introduced here by the vibrational mode (46] . Here, p(ep) stands for 
the density of states at the Fermi energy. Using the parameters of Fig. [3J, taking k — 0.1, 
we show the absolute value of these matrix elements in Fig. [5j The contour plot is mostly 
limited to values smaller than 0.1, thus we conclude that this set of parameters correspond 



to the weak coupling limit |46j. In this limit, path-integral simulations should agree with 
Master equation calculations, as we indeed confirm below. Deviations should be expected 
at larger values, k > 0.2, and we study below these cases. 

Units. We perform the simulations in arbitrary units with H = 1. One can scale all 
energies with respect to the molecule-metal hybridization T u . With T v — 1, the weak 
coupling limit covers KjY v < 0.2. To present results in physical units, we assume that all 
energy parameters are given in eV, and scale correspondingly the time unit and currents. 

Dynamics. We first focus on two representative values for the bias voltage: In the low- 
positive bias limit a stable operation is expected, reflected by a normal population, p > 
Pi. At large positive bias population inversion may take place, indicating on the onset of 
instability and potential junction rupture [271 ] . 
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FIG. 6. Population dynamics and convergence behavior of the truncated and isolated vibrational 
mode (TLS) with increasing N s . (a)-(b) Stable behavior at /zl = — zz_r = 0.2. (c)-(d) Population 
inversion at /Ul = —^jlr = 0.6. Other parameters are the same as in Fig. [3l In all figures St=l, 
N s = 3 (heavy dotted), N s = 4 (heavy dashed), N s = 5 (dashed-dotted), N s = 6 (dotted), N s = 7 
(dashed) and N s = 8 (full). We used L s = 30 electronic states at each fermionic bath with sharp 
cutoffs at ±1. 
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FIG. 7. Independence of the population po on the initial state for different biases, = 0.4, 0.8, 
1.2 top to bottom. Other parameters are the same as in Fig. [3] and Fig. [6l 



Fig. O displays the TLS dynamics, and we present data for different memory sizes N s St. 
At small positive bias, e^ — e a < ujq, the mode occupation is "normal", po > p±. In particular, 
in panels (a)-(b) we discern the case /il = — fiR = 0.2, resulting in the (shifted) electronic 
energies = and e a = 0.2. In this case the (converged) asymptotic long-time population 
(representing steady-state values), are Pq s = 0.76 and p{ s = 0.24. In contrast, when the bias 
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FIG. 8. Population dynamics, po(t). (a) Comparison between exact simulations (dashed) and 
Master equation results (dashed-dotted) at k = 0.2. (b) Deviations between exact results and 
Master equations for k = 0.1 (dot and o) and for k = 0.2 (+ and x). Other parameters are as 
determined in Fig. [3j 



is large, \xl = —fiR = 0.6, the electronic levels are shifted to e<2 = 0.4 and e a = —0.2, and 
electrons crossing the junction discard their excess energy into the vibrational mode. Indeed, 
we see in Fig. H^c)-(d) the process of population inversion, Pq s = 0.43 and pl s = 0.57. The 
TLS approaches the steady-state value around t ss ~ 0.1 ps. Regarding convergence behavior, 
we note that at large bias convergence is reached with a shorter memory size, compared to 



the small bias expected [171 ] . 



Fig. [7J exhibits the dynamics with different initial conditions, demonstrating that the 
steady-state value is identical, yet the timescale to reach the stationary limit may depend 
on the initial state. 

We compare the exact dynamics to the Master equation time evolution behavior, reached 
by solving Eq. (I52p . Panel (a) in Fig. [S] demonstrates excellent agreement for k = 0.2, 
for both positive and negative biases. Below we show that at this value Master equation's 
predictions for the charge current deviate from the exact result. Panel (b) in Fig. [8] focuses 
on the departure of Master equation data from the exact values. These deviations are small, 
but their dynamics indicate on the existence of high order excitation and relaxation rates, 
beyond the second order rates of Sec. IIVDI 

Steady-state characteristics. The full bias scan of the steady-state population is displayed 
in Fig. |9l and we compare path-integral results with Master equation calculations, revealing 
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an excellent agreement in this weak coupling limit (k = 0.1). The convergence behavior 
is presented in Fig. (TOj and we plot the steady-state values as a function of memory size 
(t c ) for three different time steps, for representative biases. The path-integral results well 
converge at intermediate-to-large positive biases, A/x > 0.2. We had difficulty converging 
our results in two domains: (i) At small-positive potential bias, A/x < 0.2. Here, large 
memory size should be used for reaching full convergence; decorrelation time approximately 
scales with 1/A/x. (ii) At large negative biases, A/x < —0.4 the current is very small as we 
show immediately. This implies poor convergence at the range of r c employed. At these 
negative biases the data oscillates with r c , thus at negative bias it is the averaged value for 
several-large r c which is plotted in Fig. 

Charge current. We show the current characteristics in Fig. (TTJ and confirm that the 
junction acts as a charge rectifier. The insets display transient data, affirming that at large 
bias steady-state is reached faster than in the low bias case. 

Strong coupling. Results at weak-to-strong couplings are shown in Fig. [12j The value 
of the current, as reached from Master equation calculations, scale with k 2 . In contrast, 
exact simulations indicate that the current grows more slowly with k, and it displays clear 
deviations (up to 50%) from the perturbative Master equation result at k — 0.3. Inter- 
estingly, the vibrational occupation (inset) shows little sensitivity to the coupling strength, 
and even at k = 0.3 the Master equation technique provides an excellent estimation for the 
levels occupation. This could be reasoned by the fact that excited levels occupation is given 
by ratio of excitation rates to the sum of excitation and relaxation rates. Such a ratio is 
(apparently) only weakly sensitive to the value of k itself, even when high-order processes 
do contribute to the current. 

Direct tunneling vs. vibrational assisted transport. Until this point (and beyond this 
subsection) we have taken Vda — 0. We now evaluate the contribution of different trans- 
port mechanisms by adding a direct D-A tunneling term, Vd a 7^ to our model Hamilto- 
nian. Electrons can now either cross the junction in a coherent manner, or inelastically, 
by exciting/de-exciting the vibrational mode. Fig. H3] demonstrates that when the vibra- 
tion assisted transport energy k is identical in strength to the direct tunneling element 
Vda, the overall current is enhanced by about a factor of two, compared to the case when 
only vibrational-assisted processes are allowed. We also note that the occupation of the 
vibrational mode is barely affected by the opening of the new electron transmission route 
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FIG. 9. Converged data for the population of the isolated vibrational mode in the steady-state 
limit with k = 0.1. Other parameters are the same as in Fig. [3l We display path-integral data for 
Po (o) and pi (□). Master equation results appear as dashed line for po and dashed-dotted line for 
Pi- 



(deviations are within the convergence error). While we compare IF data to Master equa- 
tion results when Vd a = 0, in the general case of a nonzero D-A tunneling term perturbative 
methods are more involved, and techniques similar to those developed for the AH model 



should be used 



10 



-121. 



2. Equilibration with a secondary phonon bath 

We couple the isolated-truncated vibrational mode to a secondary phonon bath, and 
follow the mode equilibration with this bath and the removal of the vibrational instability 
effect, as we increase the vibrational mode-phonon bath coupling. As an initial condition, 
the boson bath is assumed to be thermal with an inverse temperature f3 p h- This bath is 
characterized by an ohmic spectral function ( 1501) with the dimensionless Kondo parameter 
Kd, characterizing subsystem-bath coupling, and the cutoff frequency co c . 

Population behavior. We follow the mode dynamics to the steady-state limit using the 
path-integral approach of Sec. IIHI The bosonic IF is given in the appendix. We compare 
exact results with Master equation predictions, and Fig. [TJ] depicts our simulations. The 
following observations can be made: (i) The vibrational instability effect is removed already 
for Kd = 0.01, though nonequilibrium effects are still largely visible in the mode occupation, 
(ii) The vibrational mode is closed to be equilibrated with the phonon bath once Kd ~ 0.1. 
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FIG. 10. (a) Convergence behavior of the population po i n the steady-state limit for n = 0.1. 
Other parameters are the same as in Fig. [3l Plotted are the steady-state values using different 
time steps, 5t = 0.8 (o), 8t = 1.0 (□), and 5t = 1.2 (o) at different biases, as indicated at the right 
end. (b) Population mean and its standard deviation, utilizing the last six points from panel (a), 
(c) Current mean and its standard deviation, similarly attained from the data in panel (a). 
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transient results at A/x = 1.0 eV (top) and A// = —0.5 eV (bottom). 

(iii) For the present range of parameters (large u c , weak subsystem-bath couplings), Master 
equation tools reproduce the behavior of the vibrational mode. 

Charge Current. The role of the secondary phonon bath on the charge current char- 
acteristics is displayed in Fig. [13 There are two main effects related to the presence of 
the phonon bath: The step structure about zero bias is flattened when ~ 0.1, and the 
current-voltage characteristics as a whole is slightly enhanced at finite K d , at large bias. 
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FIG. 12. Charge current and vibrational occupation in the steady-state limit at different electron- 
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(□). The corresponding Master equation results appear as dashed lines. Inset: The population 
behavior in the steady-state limit for the three cases k = 0.1 (o), k = 0.2 (o) and k = 0.3 (□), with 
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FIG. 13. Study of the contribution of different transport mechanisms. Vda = (o), with Master 
equation results noted by the dashed line, and Vd a = 0.1 (□). The main plot displays the charge 
current. The inset presents the vibrational levels occupation, with empty symbols for po and filled 
symbols for p\. Other parameters are the same as in Fig. [3l particularly, the vibrational-electronic 
coupling is n = 0.1. 

Both of these effects are excellently reproduced with the Master equation, and we conclude 
that in this weak-coupling regime the presence of the phonon bath does not affect the recti- 
fying behavior of the junction. We have also verified (not shown) that at stronger coupling, 
k = 0.2 (where Master equation fails), the thermal bath similarly affects the current- voltage 
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behavior. 

An important observation is that the current itself does not testify on the state of the 
vibrational mode, whether it is in a stable or an unstable nonequilibrium state, and whether 
it is thermalized. The study of the current characteristics itself (j vs. A/i) is therefore 
insufficient to determine junction stability. More detailed information can be gained from 



the structure of the first derivative, <ij/(i(A/x), the 



derivative, d 2 j/dAfi 2 , providing spectral features [47 



ocal density of states, and the second 



49| . In order to examine these quanti- 



ties, our simulations should be performed with many more bath states, to eliminate possible 
spurious oscillations in the current (of small amplitudes) that may result from the finite 
discretization of the fermi baths. 
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FIG. 14. Equilibration of the molecular vibrational mode with increasing coupling to a secondary 
phonon bath. Path-integral results, (full symbols for pi, and empty symbols for po) with = 
(o), Kd = 0.01 (o), Kd = 0.1 (□), and, Kd = 0.1, k = (<). Unless otherwise specified, k = 0.1, 
Pph = 5 and the spectral function follows (|50p with u c =15. All other electronic parameters are the 
same as in Fig. [3j Master equation results appear in dotted lines. 



F. Convergence and Computational aspects 



Convergence of the path-integral method should be verified with respect to three numer- 
ical parameters: the number of states used to mimic a fermi sea, L s , the time step adopted, 
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FIG. 15. Charge current for an isolated mode,!^ = (o), and an equilibrated mode, = 0.1, 
(3 p h=5, uj c =15 (□). Other electronic parameters are given in Fig. [3j Master equation results 
appear in dashed-dotted lines. 

St, and the memory time accounted for, r c . (i) Fermi sea discretization. We have found that 
excellent convergence is achieved for relatively "small" fermi reservoirs, taking into account 
L s > 20 states for each reservoir. In our simulations we practically adopted L s = 30 for 
each Fermi bath, (ii) Time-step discretization. The first criteria in selecting the value of the 
time step St is that dynamical features of the isolated vibrational mode should be observed. 
Using uq = 0.2, the period of the bath-free Rabbi oscillation is 2ti/(ujo) ~ 30, thus a time 
step of St ~ 1 can capture the details of the TLS oscillation. This consideration serves as an 
"upper bound" criteria. The second consideration connects to the time discretization error, 
originates from the approximate splitting of the total time evolution operator into a product 
of terms, see Eq. ([5l). For the particular Trotter dec omp osition employed, the leading error 



grows with St 3 x ([H s , [V, H s ]]/12 + [V, [V, H s ]}/24) [50j where V = V SB + V SF + H B + H F . 



The decomposition is exact when the coupling of the subsystem to the reservoirs is weak 
and the time-step is small, St — > 0. For large coupling one should take a sufficiently small 
time-step in order to avoid significant error buildup. In the preset work, the dimensionless 
coupling to the fermi sea npnX^Xr is typically maintained lower that 0.3; the dimensionless 
coupling to the boson bath is taken as = 0.1. The value of St = 0.6 — 1.2 is thus suf- 
ficiently small for our simulations, (iii) Memory error. Our approach assumes that bath 
correlations exponentially decay resulting from the finite temperature and the nonequilib- 
rium condition. Based on this assumption, the total influence functional was truncated to 
include only a finite number of time steps N s , where r c = N s St. The total IF is retrieved 
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by taking the limit N s — > N, (N = t/5t). Our simulations were performed for N s = 3. ..9, 
covering memory time up to r c = N s 5t ~ 10. The results displayed converged for N s ~ 7 — 9 
for 5t = 1. 

Computational efforts can be partitioned into two parts: In the initialization step the 
(time invariant) IFs are computed. The size of the fermionic IF is d 2Ns , where d is the 
dimensionality of the subsystem (two in our simulations). The power of two at the exponent 
results from the forward and backward time evolution operators in the path-integral expres- 
sion. This initialization effort thus scales exponentially with the memory size accounted for. 
The preparation of the bosonic IF is more efficient if the FV IF is used [16j ]. In the second, 
time evolution, stage, we iteratively apply the linear map ( |27|) or (134"]) . a multiplication of 
two objects of length d 2Ns . This operation linearly scales in time. 

We now comment on the simulation time of a convergence analysis as presented in Fig. 
[TOl covering three different time steps and N s = 3, ...,9. The MATLAB implementation of 
the computational algorithm took advantage of the MATLAB built-in multi-threaded par- 
allel features and utilized 100% of all available CPU cores on a node. When executed on 
one cluster node with two quad-core 2.2GHz AMD Opteron cpus and 16GB memory, con- 
vergence analysis for of the full voltage scan took about 4x24 hours and 250MB of memory. 



Computations performed on the GPC supercomputer at the SciNet HPC Consortium [51] 
were three times faster. Computational time scales linearly with the simulated time t. For 
a fixed N s value, the computational effort does not depend on the system temperature and 
other parameters employed. 



V. SUMMARY 

We have developed an iterative numerically-exact path-integral scheme that can follow 
the dynamics, to the steady-state limit, of subsystems coupled to multiple bosonic and 
fermionic reservoirs in an out-of-equilibrium initial state. The method is based on the 
truncation of time correlations in the influence functional, beyond the memory time dictated 



by temperature and chemical biases. It combines two techniques: the QUAPI method 16], 



or treating the dynamics of subsystems coupled to harmonic baths, and the INFPI approach 



17] , useful for following the evolution of a subsystem when interacting with fermionic baths. 
The method is stable, efficient, and flexible, and it allows one to achieve transient and 
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steady-state data for both the reduced density matrix of the subsystem and expectation 
values of operators, such as the charge current and energy current. The method can be 
viewed as an extension of QUAPI, to incorporate fermions in the dynamics. It could be 
further expanded to include time-dependent Hamiltonians, e.g., pulsed fields. 

To demonstrate the method usability in the field of molecular conduction, we have applied 
the general scheme, and studied vibrational dynamics in a molecular rectifier setup, where 
vibrational equilibration with an additional phonon bath is allowed. Our main conclusions 
in this study are the following: (i) The vibrational instability effect disappears once the 
vibrational mode is weakly coupled (iQ ~ 0.01) to an additional phonon bath that can 
dissipate the excess energy, (ii) When ~ 0.1, the vibrational mode is equilibrated with 
the secondary phonon bath, (iii) The charge current does not testify on vibrational heating 
and instability. While we have performed those simulations using a truncated vibrational 
mode, a TLS, representing an anharmonic mode, we argue that the main characteristics of 
;he vibrational instability effect remain intact when the selected mode is made harmonic 



27]. 



Our simulations indicate that Master equation methods can excellently reproduce exact 
results at weak coupling, in the markovian limit. More significantly, Master equation tools 
can be used beyond the weak coupling limit (k ~ 0.3), if only a qualitative understanding 
of the junction behavior is inquired. One should note that our Master equation technique 
treats the D and A coupling to the metals exactly. It is perturbative only in the interaction 
of the vibrational mode to the electrons, and to other phonon degrees of freedom. In the 
case where tunneling transmission competes with phonon-assisted transport, only path- 
integral simulations were provided, as more involved Master equation methodologies should 
be developed in this case. 

Our future objectives are twofold: (i) to improve the time-evolution algorithm, and (ii) 
to employ the method for the study of other problems in molecular electronics and phonon- 
ics. By improving the methodology, we would like to extend the usability of our method to 



difficult parameter re gim es 



the memory function 
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strong coupling), e.g., by developing an equation-of- motion for 



53J . This will also allow us to simulate more feasibly the dynamics 



of an n— level subsystem. Another related objective is the study of heat current characteris- 



tics in the spin-boson molecular junction 54|. The single-bath spin-boson model displays a 



rich dynamics with a complex phase diagram. Similarly, we expect that the nonequilibrium 

35 



version, with two harmonic baths of different temperatures coupled to the TLS, will show 
complex behavior for its heat current- temperature characteristics. Recent results, obtained 
using an extension of the noninteracting blip approximation to the nonequilibrium regime 



55] , demonstrate rich behavior. Other problems that could be directed with our method 



include plexcitonics systems, as the coupling between surface plasmons and molecular exci- 
tons should be treated beyond the perturbative regime |56j. Finally, we have discussed the 
calculation of reduced density matrix and currents in the path-integral framework. It is of 
interest to generalize these expressions and gain higher order cumulants, for the study of 
current, noise, and fluctuation relations in many-body out-of-equilibrium systems. 
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APPENDIX A: TIME-DISCRETE FEYNMAN- VERNON INFLUENCE FUNC- 
TIONAL 

With the discretization of the path, the influence functional takes the form (j!7p . The 
coefficients tj^i were given in [16| and we include them here for the completeness of our 
presentation. The expressions are given here for the case of a single boson bath with the 
initial temperature l//3 p h and the spectral function J p h(u) = ^J^pi^pY^i 00 ~ u p)i Jph(y>) = 
Jph(-u), 
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„„, = I r a, W tf W)e -^), o <*-<*< iv 

7T J.^ w 2 sinh(/3 ph a;/2) 

= r I" <tJ-^^f^ (1 - e-«) , < k < N 
2tc J_ 00 w 2 smh(/3 p/l uy2) 

ifc,o = - r ^ ^^^^ sin(u;ft/4) sm^^K^-^, < k < N 
it J_ QO u 2 smh(/3 ph uj/2) 

r, NW = - du Jph ^ ] ex P^/ 2 ) s in( W 5t/4) sin(a;^/2)e-- < A;' < iV 
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